function [pval,GamBoot] = IPCA_char_pvalues(X,W,Nt,Gamma,F,ell,NBoot)
% Testing for significance of ell'th characteristic

% X is L*T
% W is L*L*T
% Gamma is L*K
% F is K*T

df = 5;
L = size(X,1);
[K, T] = size(F);
GammaRestricted = Gamma; GammaRestricted(ell,:) = 0;

d = NaN(L,T);
FitBeta = NaN(L,T);
FitBetaRestricted = NaN(L,T);
for t = 1:T-1
    Wt = W(:,:,t);
    Ft1 = F(:,t+1);
    Xt1 = X(:,t+1);
    FitBeta(:,t+1) = Wt*Gamma*Ft1;
    FitBetaRestricted(:,t+1) = Wt*GammaRestricted*Ft1;
    d(:,t+1) = Xt1 - FitBeta(:,t+1);
end

Gamma_Old = GammaRestricted; F_Old = F;
WaldBoot = NaN(NBoot,1);
GamBoot = NaN(NBoot,size(Gamma,2));

parfor b = 1:NBoot
    TimeIndex = randsample(2:T,T-1,true)';
    dBoot = [NaN(L,1) d(:,TimeIndex)];
    
    Student = sqrt((df-2)/df)*trnd(df,[1 T]); % unit variance and df degrees of freedom
    dBoot = bsxfun(@times, dBoot , Student);
    
    Xboot = FitBetaRestricted + dBoot;
    
    GammaBoot = IPCA_internal(Xboot,W,Nt,Gamma_Old,F_Old);
    GammaBoot_ell = GammaBoot(ell,:)';
    GamBoot(b,:) = GammaBoot_ell;
    WaldBoot(b) = GammaBoot_ell'*GammaBoot_ell;
end

Gamma_ell = Gamma(ell,:)';
Wald = Gamma_ell'*Gamma_ell;
pval = sum(WaldBoot>Wald)/NBoot;
return

%--------------------------------------------------------------------------
function Gamma_New = IPCA_internal(Xboot,W,Nt,Gamma_Old,F_Old)
MaxIterations = 1000;
Tolerance     = 1e-6;
tol           = 1;
iter          = 0;
while iter<=MaxIterations && tol>Tolerance
    [Gamma_New, F_New] = IPCA_base(Gamma_Old,Xboot,W,Nt);
    tol       = max([ abs(Gamma_New(:)-Gamma_Old(:)) ; abs(F_New(:)-F_Old(:)) ]);
    F_Old     = F_New;
    Gamma_Old = Gamma_New;
    iter      = iter+1;
end
return